#work out whether SNPs dn or ds

#don't use this function, crashes R
codingSNPref <- function(AllSNPref,dir.embl="/Paterson/Datafiles/grouse/embl_exon",seq.name=NULL,file.name=NULL){
	#AllSNPref is an object (dataframe) from baseCoverage
	#expects embl files in dir.embl, perhaps created by add.CDS
	
	if(is.null(seq.name)) AllSNPref[1,"cDNA"]
	if(is.null(file.name)) file.name=paste(dir.embl,"/",seq.name,"_exon.embl",sep="")
	require(GeneR)
	
	#for some reason won't read straight into buffer
	placeString(strReadEmbl(file.name,index="id"),seqno=0)

	#extract coding seqs
	grep.cds <- system(paste("grep \'CDS\'",file.name,"|sed \'s/FT[ \t]*CDS[ \t]*//\'"),intern=TRUE)
	grep.cds <- sub("join\\(","",grep.cds)
	grep.cds <- sub(")","",grep.cds)

	AllSNPref$cds <- FALSE;tmp.base$dn <- FALSE
	AllSNPref[,c('cds','dn')] <- t(sapply(X=AllSNPref[,"id"],FUN= catagoriseSNPs,grep.cds=grep.cds))
	AllSNPref
	}



#seq.name <- "cDNA_260-1"
#dir.embl <- "/Paterson/Datafiles/grouse/embl_exon"

#read sequence into buffer




#sequence 0 for original seq, seq 1 for mutated seq, seq 2 for temporary seq
translateCDS <- function(cds.str,seqno=0){
	#expects cds.str in for 1..10,20..30 etc
	#and for a sequence in buffer seqno
	tmp.from <- as.numeric(sapply(X=strsplit(cds.str,","),FUN=function(X){sub('\\.\\.[0-9]*$','',X)}))
	tmp.to <- as.numeric(sapply(X=strsplit(cds.str,","),FUN=function(X){sub('^[0-9]*\\.\\.','',X)}))
	assemble(seqno=seqno,from=tmp.from,to=tmp.to,strand=0,destSeqno=2)
	translate(2)
	}
	
mutateSeq <- function(snp.id,cds.str){
	#expects a snp.id, ie id from baseCoverage
	cat(snp.id,"  ")
	tmp.from <- as.numeric(sapply(X=strsplit(cds.str,","),FUN=function(X){sub('\\.\\.[0-9]*$','',X)}))
	tmp.to <- as.numeric(sapply(X=strsplit(cds.str,","),FUN=function(X){sub('^[0-9]*\\.\\.','',X)}))
	tmp.id <- strsplit(snp.id,"_")[[1]]
	if(substr(tmp.id[1],2,6)=="ontig"){
		mut.start <- as.numeric(tmp.id[2])
		}else{
		mut.start <- as.numeric(tmp.id[3])
		}
	
	
	mut.cds <- FALSE #test whether in coding region
	mut.dn <- FALSE
	if(any(mut.start>tmp.from&mut.start<tmp.to)){
		mut.cds <- TRUE
		}else{
		#if the mutation isn't in the coding region don't bother
		#working out anything else
		return(list(mut.cds=mut.cds,mut.dn=mut.dn))
		}

	mut.old <- sub("-","",tmp.id[6])
	mut.new <- sub("-","",tmp.id[7])
	if(nchar(mut.old)!=nchar(mut.new)){
		#if there's an indel in the coding sequence assume non-synonymous and return
		mut.dn <- TRUE
		return(list(mut.cds=mut.cds,mut.dn=mut.dn))
		}

	#should now be carrying on only if a straight swap of nucleotide(s)
	mut.len <- nchar(mut.old)
	mut.stop <- mut.start + mut.len -1
	tmp.seq <- paste(getSeq(seqno=0,strand=0,from=1,to=mut.start-1),
		mut.new,
		getSeq(seqno=0,strand=0,from=mut.stop+1,to=0),sep="")
	placeString(tmp.seq,seqno=1)
	old.prot <- translateCDS(cds.str,seqno=0)
	new.prot <- translateCDS(cds.str,seqno=1)
	if(!identical(old.prot,new.prot)) mut.dn <- TRUE
	
	list(mut.cds=mut.cds,mut.dn=mut.dn)
	}


catagoriseSNPs <- function(snp.id,grep.cds){
	tmp.mut.list <- vector(mode="list",length=length(grep.cds))
	for(i in 1:length(grep.cds)){
		tmp.mut.list[[i]] <- mutateSeq(snp.id,grep.cds[i])
		}
	c(any(sapply(tmp.mut.list,function(X){X$mut.cds})),
		any(sapply(tmp.mut.list,function(X){X$mut.dn})))

	}




